####### SCRIPT FOR: 
####### ACTIVE AND PASSIVE: TWO WAYS PARTY SYSTEMS INFLUENCE ELECTORAL OUTCOMES

rm(list=ls())
gc()
setwd("C:/Users/miros/Desktop/Research/Active and Passive - Two Ways Party Systems Influence Electoral Outcomes")

library(readxl)
#Numer of options dataset
num_opt_data <- read_excel("01_number_of_options_dataset.xlsx")
#Polarization dataset
pol_data <- read_excel("01_polarization_dataset_pair.xlsx")

data <- merge(num_opt_data, pol_data, by = "PAIRING")

rm(num_opt_data)
rm(pol_data)

#Computing absolute value for Ds1 used below
data$abs_d_s1 <- abs(data$d_s1)

#Disable scientific notations
options(scipen=999)

summary(data)



#Preparation of data files suitable for analyses
#Choice set size
data$abs_d_s1 <- abs(data$d_s1)  #computing absolute values for s1
num_opt <- data.frame(data$wgh_avg_option, data$avg_option,
                      data$PAIRING, data$CNTRY, data$abs_d_Ns, data$abs_d_s1, data$AVG_MAG,
                      data$e_democ, data$length_e_democ, data$e_migdppc, data$Alesina_Ethnic_Frac)
num_opt <- na.omit(num_opt)
#write.csv(num_opt, "cases_num_opt.csv")
summary(num_opt)

#Polarization: Comaprative Manifestos Project
POL_CMP <- data.frame(data$POL_CMP, data$wgh_avg_option, data$avg_option, data$simple, data$uppershare,
                      data$PAIRING, data$CNTRY, data$abs_d_Ns, data$abs_d_s1, data$AVG_MAG,
                      data$e_democ, data$length_e_democ, data$e_migdppc, data$Alesina_Ethnic_Frac)
POL_CMP <- na.omit(POL_CMP)
#write.csv(POL_CMP, "cases_pol_cmp.csv")
summary(POL_CMP)

#Polarization: Chapel Hill Expert Survey
POL_CHES <- data.frame(data$CHES_lrgen, data$wgh_avg_option, data$avg_option, data$simple, data$uppershare,
                      data$PAIRING, data$CNTRY, data$abs_d_Ns, data$abs_d_s1, data$AVG_MAG,
                      data$e_democ, data$length_e_democ, data$e_migdppc, data$Alesina_Ethnic_Frac)
POL_CHES <- na.omit(POL_CHES)
#write.csv(POL_CHES, "cases_pol_ches.csv")
summary(POL_CHES)

#Polarization: Party System Polarization Index
POL_DAL <- data.frame(data$POL_DALTON, data$wgh_avg_option, data$avg_option,
                       data$PAIRING, data$CNTRY, data$abs_d_Ns, data$abs_d_s1, data$AVG_MAG,
                       data$e_democ, data$length_e_democ, data$e_migdppc, data$Alesina_Ethnic_Frac)
POL_DAL <- na.omit(POL_DAL)
#write.csv(POL_DAL, "cases_pol_dalton.csv")
summary(POL_DAL)


### CHOICE SET SIZE
### Dataset: Constituency-LEvel Election Archive
library(devtools)
devtools::install_github("strengejacke/sjPlot")
install.packages("sjPlot")
library(ggplot2)
library(sjPlot)
library(sjmisc)
library(estimatr)

#Largest party seat share (s1): Choice set size
data$abs_d_s1 <- abs(data$d_s1)  #computing absolute values for s1
#Model 1: Base model for s1
opts.s1.m1 <- lm_robust(data.abs_d_s1 ~ log(data.wgh_avg_option),
                        data = num_opt, clusters = data.CNTRY, se_type = "stata")
summary(opts.s1.m1)

#Model 2: Base model + Controls for s1
opts.s1.m2 <- lm_robust(data.abs_d_s1 ~ log(data.wgh_avg_option) +
                        data.AVG_MAG +
                        data.e_democ + data.length_e_democ +
                        data.e_migdppc + data.Alesina_Ethnic_Frac, 
                        data = num_opt, clusters = data.CNTRY, se_type = "stata")
summary(opts.s1.m2)

#Figure:
opts.s1.p <- plot_model(opts.s1.m2,
                        type = "pred",
                        terms = "data.wgh_avg_option [exp]") +
  labs(title = NULL,
       x = "Weighted average number of options (ln)",
       y = expression(paste("Deviation: | ",italic(d[s[1]]), "| = log[", italic(s[1])/(italic(MS))^"-1/8", "]" ))) +
  theme(panel.grid.major = element_line(colour = "gray90"),
        panel.grid.minor = element_line(colour = "gray90"),
        panel.background = element_blank()) +
  geom_segment(aes(x = 1, xend = 300, y = 0, yend = 0)) +
  geom_segment(aes(x = 1, xend = 1, y = 0, yend = 0.5)) +
  coord_trans(limx = c(0, 35), limy = c(0, 0.301))
opts.s1.p


#Effective number of parties (Ns): Choice set size
#Model 3: Base model for Ns
opts.Ns.m1 <- lm_robust(data.abs_d_Ns ~ log(data.wgh_avg_option),
                        data = num_opt, clusters = data.CNTRY, se_type = "stata")
summary(opts.Ns.m1)

#Model 4: Base model + Controls for Ns
opts.Ns.m2 <- lm_robust(data.abs_d_Ns ~ log(data.wgh_avg_option) +
                        data.AVG_MAG +
                        data.e_democ + data.length_e_democ +
                        data.e_migdppc + data.Alesina_Ethnic_Frac, 
                        data = num_opt, clusters = data.CNTRY, se_type = "stata")
summary(opts.Ns.m2)

#Figure:
opts.Ns.p <- plot_model(opts.Ns.m2,
                             type = "pred",
                             terms = "data.wgh_avg_option [exp]") +
  labs(title = NULL,
       x = "Weighted average number of options (ln)",
       y = expression(paste("Deviation: | ",italic(d[N[s]]), "| = log[", italic(N[s])/(italic(MS))^"1/6", "]" ))) +
  theme(panel.grid.major = element_line(colour = "gray90"),
        panel.grid.minor = element_line(colour = "gray90"),
        panel.background = element_blank()) +
  geom_segment(aes(x = 1, xend = 300, y = 0, yend = 0)) +
  geom_segment(aes(x = 1, xend = 1, y = 0, yend = 0.5)) +
  coord_trans(limx = c(0, 35), limy = c(0, 0.301))
opts.Ns.p


#Multiplots
library(ggpubr)
#png(file = "Fig1_opts.png", width = 4000, height = 2000, res = 600)
ggarrange(opts.s1.p, opts.Ns.p, ncol=2, nrow=1)
#dev.off()

#Table
library(stargazer)
library(estimatr)
opts.s1.m1.lm <- lm(data.abs_d_s1 ~ log(data.wgh_avg_option), data = num_opt)
opts.s1.m2.lm <- lm(data.abs_d_s1 ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = num_opt)
opts.Ns.m1.lm <- lm(data.abs_d_Ns ~ log(data.wgh_avg_option), data = num_opt)
opts.Ns.m2.lm <- lm(data.abs_d_Ns ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = num_opt)
stargazer(opts.s1.m1.lm, opts.s1.m2.lm, opts.Ns.m1.lm, opts.Ns.m2.lm,
          se = starprep(opts.s1.m1, opts.s1.m2, opts.Ns.m1, opts.Ns.m2, 
          clusters = data$CNTRY, se_type = "stata"),
          type="html", out="model_1-4_choice-set.htm")


### POLARIZATION
#correlation among measures
library(Hmisc)
data_pol <- data.frame(data$POL_CMP, data$CHES_lrgen, data$POL_DALTON)
rcorr(as.matrix(data_pol), type="pearson")

library(devtools)
devtools::install_github("strengejacke/sjPlot")
install.packages("sjPlot")
library(ggplot2)
library(sjPlot)
library(sjmisc)
library(estimatr)

#s1: Comparative Manifesto Project
#Model 5: Base Model
pol.cmp.s1.1 <- lm_robust(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option), 
                          data = POL_CMP, clusters = data.CNTRY, se_type = "stata")
summary(pol.cmp.s1.1)

#Model 6: Base Model
pol.cmp.s1.2 <- lm_robust(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option)+
                            data.AVG_MAG +
                            data.e_democ + data.length_e_democ +
                            data.e_migdppc + data.Alesina_Ethnic_Frac,  
                          data = POL_CMP, clusters = data.CNTRY, se_type = "stata")
summary(pol.cmp.s1.2)

#Figure:
pol.cmp.s1.1.p <- plot_model(pol.cmp.s1.2,
                             type = "pred",
                             terms = "data.wgh_avg_option [exp]") +
  labs(title = "Data: The Manifesto Project",
       x = "Weighted average number of options (ln)",
       y = expression(paste("Deviation: | ",italic(d[s[1]]), "| = log[", italic(s[1])/(italic(MS))^"-1/8", "]" ))) +
  theme(panel.grid.major = element_line(colour = "gray90"),
        panel.grid.minor = element_line(colour = "gray90"),
        panel.background = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  geom_segment(aes(x = 1, xend = 300, y = 0, yend = 0)) +
  geom_segment(aes(x = 1, xend = 1, y = 0, yend = 0.5)) +
  coord_trans(limx = c(0, 35), limy = c(0, 0.301))
pol.cmp.s1.1.p


#Ns: Comparative Manifesto Project
#Model 7: Base model
pol.cmp.Ns.1 <- lm_robust(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option), 
                          data = POL_CMP, clusters = data.CNTRY, se_type = "stata")
summary(pol.cmp.Ns.1)

#Model 8: Base model + controls
pol.cmp.Ns.2 <- lm_robust(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option) +
                          data.AVG_MAG +
                          data.e_democ + data.length_e_democ +
                          data.e_migdppc + data.Alesina_Ethnic_Frac, 
                          data = POL_CMP, clusters = data.CNTRY, se_type = "stata")
summary(pol.cmp.Ns.2)

#Figure:
pol.cmp.Ns.1.p <- plot_model(pol.cmp.Ns.2,
                             type = "pred",
                             terms = "data.wgh_avg_option [exp]") +
  labs(title = "", #"Data: The Manifesto Project",
       x = "Weighted average number of options (ln)",
       y = expression(paste("Deviation: | ",italic(d[N[s]]), "| = log[", italic(N[s])/(italic(MS))^"1/6", "]" ))) +
  theme(panel.grid.major = element_line(colour = "gray90"),
        panel.grid.minor = element_line(colour = "gray90"),
        panel.background = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  geom_segment(aes(x = 1, xend = 300, y = 0, yend = 0)) +
  geom_segment(aes(x = 1, xend = 1, y = 0, yend = 0.5)) +
  coord_trans(limx = c(0, 35), limy = c(0, 0.301))
pol.cmp.Ns.1.p



#s1: Chapel Hill Expert Survey
#Model 9:
pol.ch.s1.1 <- lm_robust(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option), 
                         data = POL_CHES, clusters = data.CNTRY, se_type = "stata")
summary(pol.ch.s1.1)

#Model 10:
pol.ch.s1.2 <- lm_robust(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option) +
                         data.AVG_MAG +
                         data.e_democ + data.length_e_democ +
                         data.e_migdppc + data.Alesina_Ethnic_Frac,   
                         data = POL_CHES, clusters = data.CNTRY, se_type = "stata")
summary(pol.ch.s1.2)

pol.ch.s1.1.p <- plot_model(pol.ch.s1.2,
                            type = "pred",
                            terms = "data.wgh_avg_option [exp]") +
  labs(title = "Data: Chapel Hill Expert Survey",
       x = "Weighted average number of options (ln)",
       y = expression(paste("Deviation: | ",italic(d[s[1]]), "| = log[", italic(s[1])/(italic(MS))^"-1/8", "]" ))) +
  theme(panel.grid.major = element_line(colour = "gray90"),
        panel.grid.minor = element_line(colour = "gray90"),
        panel.background = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  geom_segment(aes(x = 1, xend = 300, y = 0, yend = 0)) +
  geom_segment(aes(x = 1, xend = 1, y = 0, yend = 0.5)) +
  coord_trans(limx = c(0, 35), limy = c(0, 0.301))
pol.ch.s1.1.p


#Ns: Chapel Hill Expert Survey
#Model 11:
pol.ch.Ns.1 <- lm_robust(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option), 
                         data = POL_CHES, clusters = data.CNTRY, se_type = "stata")
summary(pol.ch.Ns.1)

#Model 12:
pol.ch.Ns.2 <- lm_robust(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option) +
                         data.AVG_MAG +
                         data.e_democ + data.length_e_democ +
                         data.e_migdppc + data.Alesina_Ethnic_Frac,  
                         data = POL_CHES, clusters = data.CNTRY, se_type = "stata")
summary(pol.ch.Ns.2)

#Figure:
pol.ch.Ns.1.p <- plot_model(pol.ch.Ns.2,
                             type = "pred",
                             terms = "data.wgh_avg_option [exp]") +
  labs(title = "", #"Data: Chapel Hill Expert Survey",
       x = "Weighted average number of options (ln)",
       y = expression(paste("Deviation: | ",italic(d[N[s]]), "| = log[", italic(N[s])/(italic(MS))^"1/6", "]" ))) +
  theme(panel.grid.major = element_line(colour = "gray90"),
        panel.grid.minor = element_line(colour = "gray90"),
        panel.background = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  geom_segment(aes(x = 1, xend = 300, y = 0, yend = 0)) +
  geom_segment(aes(x = 1, xend = 1, y = 0, yend = 0.5)) +
  coord_trans(limx = c(0, 35), limy = c(0, 0.301))
pol.ch.Ns.1.p


#s1: Party System Polarization Index
#Model 13:
pol.dal.s1.1 <-  lm_robust(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option), 
                           data = subset(POL_DAL, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
summary(pol.dal.s1.1)

#Model 14:
pol.dal.s1.2 <-  lm_robust(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option)+
                           data.AVG_MAG +
                           data.e_democ + data.length_e_democ +
                           data.e_migdppc + data.Alesina_Ethnic_Frac,  
                           data = subset(POL_DAL, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
summary(pol.dal.s1.2)

#Figure:
pol.dal.s1.1.p <- plot_model(pol.dal.s1.2,
                             type = "pred",
                             terms = "data.wgh_avg_option [exp]") +
  labs(title = "Data: Party System Polarization Index",
       x = "Weighted average number of options (ln)",
       y = expression(paste("Deviation: | ",italic(d[s[1]]), "| = log[", italic(s[1])/(italic(MS))^"-1/8", "]" ))) +
  theme(panel.grid.major = element_line(colour = "gray90"),
        panel.grid.minor = element_line(colour = "gray90"),
        panel.background = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  geom_segment(aes(x = 1, xend = 300, y = 0, yend = 0)) +
  geom_segment(aes(x = 1, xend = 1, y = 0, yend = 0.5)) +
  coord_trans(limx = c(0, 35), limy = c(0, 0.301))
pol.dal.s1.1.p


#Ns: Party Polarization Index
#Model 15:
pol.dal.Ns.1 <- lm_robust(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option), 
                          data = subset(POL_DAL, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
summary(pol.dal.Ns.1)

#Model 16:
pol.dal.Ns.2 <- lm_robust(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option) +
                          data.AVG_MAG +
                          data.e_democ + data.length_e_democ +
                          data.e_migdppc + data.Alesina_Ethnic_Frac,  
                          data = subset(POL_DAL, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
summary(pol.dal.Ns.2)

#Figure:
pol.dal.Ns.1.p <- plot_model(pol.dal.Ns.2,
                            type = "pred",
                            terms = "data.wgh_avg_option [exp]") +
  labs(title = "", #"Data: Party System Polarization Index",
       x = "Weighted average number of options (ln)",
       y = expression(paste("Deviation: | ",italic(d[N[s]]), "| = log[", italic(N[s])/(italic(MS))^"1/6", "]" ))) +
  theme(panel.grid.major = element_line(colour = "gray90"),
        panel.grid.minor = element_line(colour = "gray90"),
        panel.background = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  geom_segment(aes(x = 1, xend = 300, y = 0, yend = 0)) +
  geom_segment(aes(x = 1, xend = 1, y = 0, yend = 0.5)) +
  coord_trans(limx = c(0, 35), limy = c(0, 0.301))
pol.dal.Ns.1.p



#Multiplots
library(ggpubr)
#png(file = "Fig2_pol.png", width = 7000, height = 4000, res = 600)
ggarrange(pol.cmp.s1.1.p, pol.ch.s1.1.p, pol.dal.s1.1.p, 
          pol.cmp.Ns.1.p, pol.ch.Ns.1.p, pol.dal.Ns.1.p,
          ncol=3, nrow=2)
#dev.off()

#Tables
library(stargazer)
library(estimatr)

#Table: Manifestos project
#Running lm() models for coefficients
pol.cmp.s1.1.lm <- lm(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option), data = POL_CMP)
pol.cmp.s1.2.lm <- lm(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP)
pol.cmp.Ns.1.lm <- lm(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option), data = POL_CMP)
pol.cmp.Ns.2.lm <- lm(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP)
#Table:
stargazer(pol.cmp.s1.1.lm, pol.cmp.s1.2.lm, pol.cmp.Ns.1.lm, pol.cmp.Ns.2.lm,
          se = starprep(pol.cmp.s1.1, pol.cmp.s1.2, pol.cmp.Ns.1, pol.cmp.Ns.2, 
          clusters = POL_CMP$CNTRY, se_type = "stata"),
          type="html", out="model_5-8_pol_cmp.htm")

#Table: Chapel Hill Expert Survey
#Running lm() models for coefficients
pol.ch.s1.1.lm <- lm(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option), data = POL_CHES)
pol.ch.s1.2.lm <- lm(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CHES)
pol.ch.Ns.1.lm <- lm(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option), data = POL_CHES)
pol.ch.Ns.2.lm <- lm(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CHES)
#Table:
stargazer(pol.ch.s1.1.lm, pol.ch.s1.2.lm, pol.ch.Ns.1.lm, pol.ch.Ns.2.lm,
          se = starprep(pol.ch.s1.1, pol.ch.s1.2, pol.ch.Ns.1, pol.ch.Ns.2, 
                        clusters = POL_CHES$CNTRY, se_type = "stata"),
          type="html", out="model_9-12_pol_ches.htm")

#Table: Party System Polarization Index
#Running lm() models for coefficients
pol.dal.s1.1.lm <- lm(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option), data = subset(POL_DAL, data.POL_DALTON != 0))
pol.dal.s1.2.lm <- lm(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = subset(POL_DAL, data.POL_DALTON != 0))
pol.dal.Ns.1.lm <- lm(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option), data = subset(POL_DAL, data.POL_DALTON != 0))
pol.dal.Ns.2.lm <- lm(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = subset(POL_DAL, data.POL_DALTON != 0))
#Table:
stargazer(pol.dal.s1.1.lm, pol.dal.s1.2.lm, pol.dal.Ns.1.lm, pol.dal.Ns.2.lm,
          se = starprep(pol.dal.s1.1, pol.dal.s1.2, pol.dal.Ns.1, pol.dal.Ns.2,
                        clusters = POL_DAL$CNTRY, se_type = "stata"),
          type="html", out="model_13-16_pol_dal.htm")



##########
# APPENDIX: Robustness checks

#Preparation of data file
#Choice set size
data$abs_d_s1 <- abs(data$d_s1)  #computing absolute values for s1
num_opt_rc <- data.frame(data$wgh_avg_option, data$AVG_MAG,
                         data$PAIRING, data$CNTRY, data$abs_d_Ns, data$abs_d_s1, 
                         data$e_democ, data$length_e_democ, data$e_fh_ipolity2,
                         data$e_migdppc,
                         data$Alesina_Ethnic_Frac, data$Fearon_Ethnic_Frac)
num_opt_rc <- na.omit(num_opt_rc)
#write.csv(num_opt_rc, "cases_num_opt.csv")

#s1 vs choice set size
library(estimatr)
num_rc.s1.m1 <- lm_robust(data.abs_d_s1 ~ log(data.wgh_avg_option), data = num_opt_rc, clusters = data.CNTRY, se_type = "stata")
num_rc.s1.m2 <- lm_robust(data.abs_d_s1 ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = num_opt_rc, clusters = data.CNTRY, se_type = "stata")
num_rc.s1.m3 <- lm_robust(data.abs_d_s1 ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = num_opt_rc, clusters = data.CNTRY, se_type = "stata")
num_rc.s1.m4 <- lm_robust(data.abs_d_s1 ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = num_opt_rc, clusters = data.CNTRY, se_type = "stata")
num_rc.s1.m5 <- lm_robust(data.abs_d_s1 ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = num_opt_rc, clusters = data.CNTRY, se_type = "stata")
num_rc.s1.m1.lm <- lm(data.abs_d_s1 ~ log(data.wgh_avg_option), data = num_opt_rc)
num_rc.s1.m2.lm <- lm(data.abs_d_s1 ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = num_opt_rc)
num_rc.s1.m3.lm <- lm(data.abs_d_s1 ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = num_opt_rc)
num_rc.s1.m4.lm <- lm(data.abs_d_s1 ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = num_opt_rc)
num_rc.s1.m5.lm <- lm(data.abs_d_s1 ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = num_opt_rc)

library(stargazer)
stargazer(num_rc.s1.m1.lm, num_rc.s1.m2.lm, num_rc.s1.m3.lm, num_rc.s1.m4.lm, num_rc.s1.m5.lm,
          se = starprep(num_rc.s1.m1, num_rc.s1.m2, num_rc.s1.m3, num_rc.s1.m4, num_rc.s1.m5,
                        clusters = num_opt_rc$data.CNTRY, se_type = "stata"),
          type="html", out="rc01_num-opt_s1.htm")

#Ns vs choice set size
library(estimatr)
num_rc.Ns.m1 <- lm_robust(data.abs_d_Ns ~ log(data.wgh_avg_option), data = num_opt_rc, clusters = data.CNTRY, se_type = "stata")
num_rc.Ns.m2 <- lm_robust(data.abs_d_Ns ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = num_opt_rc, clusters = data.CNTRY, se_type = "stata")
num_rc.Ns.m3 <- lm_robust(data.abs_d_Ns ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = num_opt_rc, clusters = data.CNTRY, se_type = "stata")
num_rc.Ns.m4 <- lm_robust(data.abs_d_Ns ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = num_opt_rc, clusters = data.CNTRY, se_type = "stata")
num_rc.Ns.m5 <- lm_robust(data.abs_d_Ns ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = num_opt_rc, clusters = data.CNTRY, se_type = "stata")
num_rc.Ns.m1.lm <- lm(data.abs_d_Ns ~ log(data.wgh_avg_option), data = num_opt_rc)
num_rc.Ns.m2.lm <- lm(data.abs_d_Ns ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = num_opt_rc)
num_rc.Ns.m3.lm <- lm(data.abs_d_Ns ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = num_opt_rc)
num_rc.Ns.m4.lm <- lm(data.abs_d_Ns ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = num_opt_rc)
num_rc.Ns.m5.lm <- lm(data.abs_d_Ns ~ log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = num_opt_rc)

library(stargazer)
stargazer(num_rc.Ns.m1.lm, num_rc.Ns.m2.lm, num_rc.Ns.m3.lm, num_rc.Ns.m4.lm, num_rc.Ns.m5.lm,
          se = starprep(num_rc.Ns.m1, num_rc.Ns.m2, num_rc.Ns.m3, num_rc.Ns.m4, num_rc.Ns.m5,
                        clusters = num_opt_rc$data.CNTRY, se_type = "stata"),
          type="html", out="rc02_num-opt_Ns.htm")


#Polarization: Manifesto Project
POL_CMP_rc <- data.frame(data$POL_CMP, data$wgh_avg_option, data$AVG_MAG,
                      data$PAIRING, data$CNTRY, data$abs_d_Ns, data$abs_d_s1, 
                      data$e_democ, data$length_e_democ, data$e_fh_ipolity2,
                      data$e_migdppc,
                      data$Alesina_Ethnic_Frac, data$Fearon_Ethnic_Frac)
POL_CMP_rc <- na.omit(POL_CMP_rc)
#write.csv(POL_CMP_rc, "cases_pol_cmp.csv")

#s1 vs Manifesto Project
library(estimatr)
pol_cmp.s1.m1 <- lm_robust(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option), data = POL_CMP_rc, clusters = data.CNTRY, se_type = "stata")
pol_cmp.s1.m2 <- lm_robust(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP_rc, clusters = data.CNTRY, se_type = "stata")
pol_cmp.s1.m3 <- lm_robust(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CMP_rc, clusters = data.CNTRY, se_type = "stata")
pol_cmp.s1.m4 <- lm_robust(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP_rc, clusters = data.CNTRY, se_type = "stata")
pol_cmp.s1.m5 <- lm_robust(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CMP_rc, clusters = data.CNTRY, se_type = "stata")
pol_cmp.s1.m1.lm <- lm(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option), data = POL_CMP_rc)
pol_cmp.s1.m2.lm <- lm(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP_rc)
pol_cmp.s1.m3.lm <- lm(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CMP_rc)
pol_cmp.s1.m4.lm <- lm(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP_rc)
pol_cmp.s1.m5.lm <- lm(data.abs_d_s1 ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CMP_rc)

library(stargazer)
stargazer(pol_cmp.s1.m1.lm, pol_cmp.s1.m2.lm, pol_cmp.s1.m3.lm, pol_cmp.s1.m4.lm, pol_cmp.s1.m5.lm,
          se = starprep(pol_cmp.s1.m1, pol_cmp.s1.m2, pol_cmp.s1.m3, pol_cmp.s1.m4, pol_cmp.s1.m5,
                        clusters = POL_CMP_rc$data.CNTRY, se_type = "stata"),
          type="html", out="rc03_pol-cmp_s1.htm")

#Ns vs Manifesto Project
library(estimatr)
pol_cmp.Ns.m1 <- lm_robust(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option), data = POL_CMP_rc, clusters = data.CNTRY, se_type = "stata")
pol_cmp.Ns.m2 <- lm_robust(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP_rc, clusters = data.CNTRY, se_type = "stata")
pol_cmp.Ns.m3 <- lm_robust(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CMP_rc, clusters = data.CNTRY, se_type = "stata")
pol_cmp.Ns.m4 <- lm_robust(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP_rc, clusters = data.CNTRY, se_type = "stata")
pol_cmp.Ns.m5 <- lm_robust(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CMP_rc, clusters = data.CNTRY, se_type = "stata")
pol_cmp.Ns.m1.lm <- lm(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option), data = POL_CMP_rc)
pol_cmp.Ns.m2.lm <- lm(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP_rc)
pol_cmp.Ns.m3.lm <- lm(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CMP_rc)
pol_cmp.Ns.m4.lm <- lm(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP_rc)
pol_cmp.Ns.m5.lm <- lm(data.abs_d_Ns ~ data.POL_CMP + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CMP_rc)

library(stargazer)
stargazer(pol_cmp.Ns.m1.lm, pol_cmp.Ns.m2.lm, pol_cmp.Ns.m3.lm, pol_cmp.Ns.m4.lm, pol_cmp.Ns.m5.lm,
          se = starprep(pol_cmp.Ns.m1, pol_cmp.Ns.m2, pol_cmp.Ns.m3, pol_cmp.Ns.m4, pol_cmp.Ns.m5,
                        clusters = POL_CMP_rc$data.CNTRY, se_type = "stata"),
          type="html", out="rc04_pol-cmp_Ns.htm")


#Polarization: Chapel Hill Expert Survey
POL_CHES_rc <- data.frame(data$CHES_lrgen, data$wgh_avg_option, data$AVG_MAG,
                       data$PAIRING, data$CNTRY, data$abs_d_Ns, data$abs_d_s1, 
                       data$e_democ, data$length_e_democ, data$e_fh_ipolity2,
                       data$e_migdppc,
                       data$Alesina_Ethnic_Frac, data$Fearon_Ethnic_Frac)
POL_CHES_rc <- na.omit(POL_CHES_rc)
#write.csv(POL_CHES, "cases_pol_ches.csv")

#s1 vs Manifesto Project
library(estimatr)
pol_ches.s1.m1 <- lm_robust(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option), data = POL_CHES_rc, clusters = data.CNTRY, se_type = "stata")
pol_ches.s1.m2 <- lm_robust(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CHES_rc, clusters = data.CNTRY, se_type = "stata")
pol_ches.s1.m3 <- lm_robust(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CHES_rc, clusters = data.CNTRY, se_type = "stata")
pol_ches.s1.m4 <- lm_robust(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CHES_rc, clusters = data.CNTRY, se_type = "stata")
pol_ches.s1.m5 <- lm_robust(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CHES_rc, clusters = data.CNTRY, se_type = "stata")
pol_ches.s1.m1.lm <- lm(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option), data = POL_CHES_rc)
pol_ches.s1.m2.lm <- lm(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CHES_rc)
pol_ches.s1.m3.lm <- lm(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CHES_rc)
pol_ches.s1.m4.lm <- lm(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CHES_rc)
pol_ches.s1.m5.lm <- lm(data.abs_d_s1 ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CHES_rc)

library(stargazer)
stargazer(pol_ches.s1.m1.lm, pol_ches.s1.m2.lm, pol_ches.s1.m3.lm, pol_ches.s1.m4.lm, pol_ches.s1.m5.lm,
          se = starprep(pol_ches.s1.m1, pol_ches.s1.m2, pol_ches.s1.m3, pol_ches.s1.m4, pol_ches.s1.m5,
                        clusters = POL_CHES_rc$data.CNTRY, se_type = "stata"),
          type="html", out="rc05_pol-ches_s1.htm")

#Ns vs Manifesto Project
library(estimatr)
pol_ches.Ns.m1 <- lm_robust(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option), data = POL_CHES_rc, clusters = data.CNTRY, se_type = "stata")
pol_ches.Ns.m2 <- lm_robust(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CHES_rc, clusters = data.CNTRY, se_type = "stata")
pol_ches.Ns.m3 <- lm_robust(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CHES_rc, clusters = data.CNTRY, se_type = "stata")
pol_ches.Ns.m4 <- lm_robust(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CHES_rc, clusters = data.CNTRY, se_type = "stata")
pol_ches.Ns.m5 <- lm_robust(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CHES_rc, clusters = data.CNTRY, se_type = "stata")
pol_ches.Ns.m1.lm <- lm(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option), data = POL_CHES_rc)
pol_ches.Ns.m2.lm <- lm(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CHES_rc)
pol_ches.Ns.m3.lm <- lm(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CHES_rc)
pol_ches.Ns.m4.lm <- lm(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CHES_rc)
pol_ches.Ns.m5.lm <- lm(data.abs_d_Ns ~ data.CHES_lrgen + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = POL_CHES_rc)

library(stargazer)
stargazer(pol_ches.Ns.m1.lm, pol_ches.Ns.m2.lm, pol_ches.Ns.m3.lm, pol_ches.Ns.m4.lm, pol_ches.Ns.m5.lm,
          se = starprep(pol_ches.Ns.m1, pol_ches.Ns.m2, pol_ches.Ns.m3, pol_ches.Ns.m4, pol_ches.Ns.m5,
                        clusters = POL_CHES_rc$data.CNTRY, se_type = "stata"),
          type="html", out="rc06_pol-ches_Ns.htm")



#Polarization: Party System Polarization Index
POL_DAL_rc <- data.frame(data$POL_DALTON, data$wgh_avg_option, data$AVG_MAG,
                      data$PAIRING, data$CNTRY, data$abs_d_Ns, data$abs_d_s1, 
                      data$e_democ, data$length_e_democ, data$e_fh_ipolity2,
                      data$e_migdppc,
                      data$Alesina_Ethnic_Frac, data$Fearon_Ethnic_Frac)
POL_DAL_rc <- na.omit(POL_DAL_rc)
#write.csv(POL_DAL, "cases_pol_dalton.csv")

#s1 vs Party System Polarization Index
library(estimatr)
pol_dal.s1.m1 <- lm_robust(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option), data = subset(POL_DAL_rc, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
pol_dal.s1.m2 <- lm_robust(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
pol_dal.s1.m3 <- lm_robust(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
pol_dal.s1.m4 <- lm_robust(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
pol_dal.s1.m5 <- lm_robust(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
pol_dal.s1.m1.lm <- lm(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option), data = subset(POL_DAL_rc, data.POL_DALTON != 0))
pol_dal.s1.m2.lm <- lm(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0))
pol_dal.s1.m3.lm <- lm(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0))
pol_dal.s1.m4.lm <- lm(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0))
pol_dal.s1.m5.lm <- lm(data.abs_d_s1 ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0))

library(stargazer)
stargazer(pol_dal.s1.m1.lm, pol_dal.s1.m2.lm, pol_dal.s1.m3.lm, pol_dal.s1.m4.lm, pol_dal.s1.m5.lm,
          se = starprep(pol_dal.s1.m1, pol_dal.s1.m2, pol_dal.s1.m3, pol_dal.s1.m4, pol_dal.s1.m5,
                        clusters = POL_DAL_rc$data.CNTRY, se_type = "stata"),
          type="html", out="rc07_pol-dal_s1.htm")

#Ns vs Party System Polarization Index
library(estimatr)
pol_dal.Ns.m1 <- lm_robust(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option), data = subset(POL_DAL_rc, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
pol_dal.Ns.m2 <- lm_robust(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
pol_dal.Ns.m3 <- lm_robust(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
pol_dal.Ns.m4 <- lm_robust(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
pol_dal.Ns.m5 <- lm_robust(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
pol_dal.Ns.m1.lm <- lm(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option), data = subset(POL_DAL_rc, data.POL_DALTON != 0))
pol_dal.Ns.m2.lm <- lm(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0))
pol_dal.Ns.m3.lm <- lm(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0))
pol_dal.Ns.m4.lm <- lm(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0))
pol_dal.Ns.m5.lm <- lm(data.abs_d_Ns ~ data.POL_DALTON + log(data.wgh_avg_option) + data.AVG_MAG + data.e_fh_ipolity2 + data.length_e_democ + data.e_migdppc + data.Fearon_Ethnic_Frac, data = subset(POL_DAL_rc, data.POL_DALTON != 0))

library(stargazer)
stargazer(pol_dal.Ns.m1.lm, pol_dal.Ns.m2.lm, pol_dal.Ns.m3.lm, pol_dal.Ns.m4.lm, pol_dal.Ns.m5.lm,
          se = starprep(pol_dal.Ns.m1, pol_dal.Ns.m2, pol_dal.Ns.m3, pol_dal.Ns.m4, pol_dal.Ns.m5,
                        clusters = POL_DAL_rc$data.CNTRY, se_type = "stata"),
          type="html", out="rc08_pol-dal_Ns.htm")


### RESPONSE TO REVIEWERS
# Link between proportionality and choice set size
ggplot(num_opt, aes(y = log(data.AVG_MAG), x = log(data.wgh_avg_option))) +
  geom_point()

rr01.mag_opt <- lm_robust(log(num_opt$data.wgh_avg_option) ~ log(num_opt$data.AVG_MAG), clusters = num_opt$data.CNTRY, se_type = "stata")
rr01.mag_opt.lm <- lm(log(num_opt$data.wgh_avg_option) ~ log(num_opt$data.AVG_MAG))

stargazer(rr01.mag_opt.lm, 
          se = starprep(rr01.mag_opt, clusters = num_opt$data.CNTRY, se_type = "stata"),
          type="html", out="rr01_mag-vs_opt.htm")

# Interaction: Proportionality vs effective number of parties
fCMP <- ggplot(POL_CMP, aes(x = data.POL_CMP, y = log(data.wgh_avg_option))) +
  geom_point() +
  geom_segment(aes(x = 0, xend = 300, y = 0, yend = 0)) +
  geom_segment(aes(x = 0, xend = 0, y = 0, yend = 7.5)) +
  coord_trans(limx = c(0, 135), limy = c(0, 5.5)) +
  labs(title = "Data: The Manifesto Project",
       y = "Weighted average number of options (ln)",
       x = "Party System Polarization\n(The Manifesto Project)") +
  theme(#plot.title = element_blank(),
        panel.grid.major.x = element_line(colour = "gray90"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(colour = "gray90"),
        panel.grid.minor.y = element_blank(),
        panel.background = element_blank())
fCMP

fCHES <- ggplot(POL_CHES, aes(x = data.CHES_lrgen, y = log(data.wgh_avg_option))) +
  geom_point() +
  geom_segment(aes(x = 0, xend = 300, y = 0, yend = 0)) +
  geom_segment(aes(x = 0, xend = 0, y = 0, yend = 7.5)) +
  coord_trans(limx = c(0, 10), limy = c(0, 5.5)) +
  labs(title = "Data: Chapel Hill Expert Survey",
       y = "Weighted average number of options (ln)",
       x = "Party System Polarization\n(Chapel Hill Expert Survey)") +
  theme(#plot.title = element_blank(),
    panel.grid.major.x = element_line(colour = "gray90"),
    panel.grid.minor.x = element_blank(),
    panel.grid.major.y = element_line(colour = "gray90"),
    panel.grid.minor.y = element_blank(),
    panel.background = element_blank())
fCHES

fDAL <- ggplot(POL_DAL, aes(x = data.POL_DALTON, y = log(data.wgh_avg_option))) +
  geom_point() +
  geom_segment(aes(x = 0, xend = 300, y = 0, yend = 0)) +
  geom_segment(aes(x = 0, xend = 0, y = 0, yend = 7.5)) +
  coord_trans(limx = c(0, 6), limy = c(0, 5.5)) +
  labs(title = "Data: Party System Polarization Index",
       y = "Weighted average number of options (ln)",
       x = "Party System Polarization\n(Party System Polarization Index)") +
  theme(#plot.title = element_blank(),
    panel.grid.major.x = element_line(colour = "gray90"),
    panel.grid.minor.x = element_blank(),
    panel.grid.major.y = element_line(colour = "gray90"),
    panel.grid.minor.y = element_blank(),
    panel.background = element_blank())
fDAL

#Multiplots
library(ggpubr)
#png(file = "FigRR_polvsfrag.png", width = 7000, height = 2300, res = 600)
ggarrange(fCMP, fCHES, fDAL,
          ncol = 3, nrow = 1)
#dev.off()

#Ns vs Party System Polarization Index
library(estimatr)
fCMP.m1 <- lm_robust(log(data.wgh_avg_option) ~ data.POL_CMP , data = POL_CMP, clusters = data.CNTRY, se_type = "stata")
fCMP.m1.lm <- lm(log(data.wgh_avg_option) ~ data.POL_CMP , data = POL_CMP)

fCHES.m1 <- lm_robust(log(data.wgh_avg_option) ~ data.CHES_lrgen , data = POL_CHES, clusters = data.CNTRY, se_type = "stata")
fCHES.m1.lm <- lm(log(data.wgh_avg_option) ~ data.CHES_lrgen , data = POL_CHES)

fDAL.m1 <- lm_robust(log(data.wgh_avg_option) ~ data.POL_DALTON , data = subset(POL_DAL_rc, data.POL_DALTON != 0), clusters = data.CNTRY, se_type = "stata")
fDAL.m1.lm <- lm(log(data.wgh_avg_option) ~ data.POL_DALTON , data = subset(POL_DAL_rc, data.POL_DALTON != 0))

library(stargazer)
stargazer(fCMP.m1.lm, fCHES.m1.lm, fDAL.m1.lm, 
          se = starprep(fCMP.m1, fCHES.m1, fDAL.m1, clusters = POL_DAL_rc$data.CNTRY, se_type = "stata"),
          type="html", out="RR_polvsfrag.htm")




#CMP
pol.cmp.s1.RR.1 <- lm_robust(data.abs_d_s1 ~ data.POL_CMP * log(data.wgh_avg_option), data = POL_CMP, clusters = data.CNTRY, se_type = "stata")
pol.cmp.s1.RR.1.lm <- lm(data.abs_d_s1 ~ data.POL_CMP * log(data.wgh_avg_option), data = POL_CMP)
pol.cmp.s1.RR.2 <- lm_robust(data.abs_d_s1 ~ data.POL_CMP * log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP, clusters = data.CNTRY, se_type = "stata")
pol.cmp.s1.RR.2.lm <- lm(data.abs_d_s1 ~ data.POL_CMP * log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP)
pol.cmp.Ns.RR.1 <- lm_robust(data.abs_d_Ns ~ data.POL_CMP * log(data.wgh_avg_option), data = POL_CMP, clusters = data.CNTRY, se_type = "stata")
pol.cmp.Ns.RR.1.lm <- lm(data.abs_d_Ns ~ data.POL_CMP * log(data.wgh_avg_option), data = POL_CMP)
pol.cmp.Ns.RR.2 <- lm_robust(data.abs_d_Ns ~ data.POL_CMP * log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP, clusters = data.CNTRY, se_type = "stata")
pol.cmp.Ns.RR.2.lm <- lm(data.abs_d_Ns ~ data.POL_CMP * log(data.wgh_avg_option) + data.AVG_MAG + data.e_democ + data.length_e_democ + data.e_migdppc + data.Alesina_Ethnic_Frac, data = POL_CMP)

library(stargazer)
stargazer(pol.cmp.s1.RR.1.lm, pol.cmp.s1.RR.2.lm, pol.cmp.Ns.RR.1.lm, pol.cmp.Ns.RR.2.lm,
          se = starprep(pol.cmp.s1.RR.1, pol.cmp.s1.RR.2, pol.cmp.Ns.RR.1, pol.cmp.Ns.RR.2, clusters = POL_CMP$data.CNTRY, se_type = "stata"),
          type="html", out="RR_CMP_polvsfrag.htm")

#Figure:
pol.cmp.s1.1.p <- plot_model(pol.cmp.s1.RR.2,
                             type = "pred",
                             terms = c("data.wgh_avg_option [exp]", "data.POL_CMP")) +
  labs(title = "Data: The Manifesto Project",
       x = "Weighted average number of options (ln)",
       y = expression(paste("Deviation: | ",italic(d[s[1]]), "| = log[", italic(s[1])/(italic(MS))^"-1/8", "]" ))) +
  theme(panel.grid.major = element_line(colour = "gray90"),
        panel.grid.minor = element_line(colour = "gray90"),
        panel.background = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  geom_segment(aes(x = 1, xend = 300, y = 0, yend = 0)) +
  geom_segment(aes(x = 1, xend = 1, y = 0, yend = 0.5)) +
  coord_trans(limx = c(0, 35), limy = c(0, 0.301))
pol.cmp.s1.1.p
